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ABSTRACT 

This paper summarizes the technical development and validation of a multiphase computational fluid 
dynamics (CFD) numerical method using the volume-of-fluid (VOF) model and a Lagrangian tracking model 
which can be employed to analyze general multiphase flow problems with free surface mechanism. The gas- 
liquid interface mass, momentum and energy conservation relationships are modeled by continuum surface 
mechanisms. A new solution method is developed such that the present VOF model can be applied for all-speed 
flow regimes. The objectives of the present study are to develop and verify the fractional volume-of-fluid cell 
partitioning approach into a predictor-corrector algorithm and to demonstrate the effectiveness of the present 
approach by simulating benchmark problems including laminar impinging jets, shear coaxial jet atomization and 
shear coaxial spray combustion flows. 


INTRODUCTION 

The atomization, breakup and combustion processes of liquid propellant in many space transportation and 
propulsion systems are of vital importance to the combustion performance and stability of the system. Numerical 
modeling of these flows poses a challenge because it requires simultaneous resolution of liquid-gas-droplets 
dynamics, and the flow regimes considered can range from very low Mach number to high-speed compressible 
flows. The flow-process modeling is further complicated by surface tension, interfacial heat and mass transfer, 
material property variation, spray formation and turbulence, and their interactions. It is the aim of this study that 
a general and robust design analysis tool will be available as a result of the incorporation of advanced numerical 
and physical models in a computational fluid dynamics flow solver - FDNS code 1 ' 2 . 

Conventionally, there have been some successful numerical algorithms for performing multiphase flow 
calculations involving free interfaces. In general, two basic approaches can be identified. One is based on the 
Lagrangian method and the other on the Eulerian frame work. The Lagrangian approach is to conform a mesh 
system with the sharp interface between gas and liquid and to track the movement of the interface according to 
sets of simplified equations governing the liquid and gas motions separately 3 . The entire Lagrangian grid 
network moves with the flee surface which makes the grid generation or remeshing procedure very complex and 
CPU intensive. In the Eulerian approach, a fixed mesh system is used to resolve the sharp interface by the 
concept of a fractional volume of fluid (VOF) 4 . The second approach is more flexible and efficient in that it 
allows the resolution of the free surface to become part of the solution. The physical processes across the sharp 
interface such as surface tension forces, droplet or wave breakup, interpbase mass transfer and 
condensation/vaporization, etc. can be modeled through the concept of a continuum surface force (CSF) method 5 
or by the method of surface reconstruction as used in the ARICC3D code 6 , SOLA-VOF code 4 . The CSF method 
is adopted in the present model for its generality and computational efficiency, especially for 3-dimensional 
computations. 

Historically, developments of numerical methods for multiphase free surface flows were primarily aimed at 
incompressible flows utilizing the Volume of Fluid (VOF) method 7, 8-9 . These methods are not effective for high- 
speed, high-Reynolds-number flows. So far, multiphase free surface flows have been considered in only a few 
general all-speed flow algorithms 6, 10, u . In the ARICC-3D Code 10 , the VOF method was implemented in the 
ALE-ICE (Arbitrary Lagrangian Eulerian-Implicit Continuous-fluid Eulerian method) algorithm for injector 
flow simulation with atomization. Due to the inefficiency of the ALE-ICE method, the ARICC-3D was found 
very time-consuming for multi-element computations. The ALE-ICE was suggested to be replaced by a pressure- 
based SIMPLE algorithm 7 which may inprove the computational efficentcy. In the RIPPLE computer program 11 , 



only finite difference solutions to the incompressible Navier-Stokes equation are obtained on an Eulerian, 
rectangular mesh using a time-split two-step projection algorithm. Compressible fluids were not solved, while 
the surface tension is treated using the CSF method. 

Based on a preliminary study of ref. 12, a unified multi-phase numerical method was developed. The gas- 
liquid interface mass, momentum and energy conservation properties are modeled by continuum surface 
mechanisms. A new solution method is developed such that the present VOF model can be applied for all-speed 
range flows. Numerical models for treating droplet dynamics are described in ref. 13. The main feature of the 
present method is to combine the advanced interface models of the Volume of Fluid (VOF) method and the 
Eulerian/Lagrangian method used in a pressure-based particulate two-phase flow solver 1 ' x 13 into a unified 
algorithm for efficient non-iterative, time-accurate calculations of multiphase free surface flows valid at all 
speeds. The present method reformulates the VOF equation to strongly couple two distinct phases (liquid and 
gas), and tracks droplets on a Lagrangian frame when a spray model is required, using a unified predictor- 
corrector technique to account for the non-linear linkages through the convective contributions of the VOF 
function. The CSF model is employed to model the surface tension forces. Formations of droplets and tracking 
of droplet dynamics are handled through the same predictor-corrector solution procedure. Thus the new 
algorithm is non-iterative and is flexible in the handling of any general geometries with arbitrarily complex 
topology in free surfaces. The present method can be applied for transient and steady state computations and can 
be implemented into any pressure-based solution method. 


GOVERNING EQUATION 


VOF methods have been mainly developed and used for low-speed flows such that incompressibility can be 
assumed. The incompressible flow assumption limits their capability. To generalize, the present formulation is 
based on compressible flow governing equations. The forms of the equations are then continuously reduced to 
their incompressible forms according to the local flow conditions and the VOF solutions. This is the uniqueness 
of the present method. To illustrate this, the density-weighted averaged conservation equation of mass, Navier- 
Stokes, and scalar variables in an Eulerian frame work can be written as: 
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where D and S represent diffusion and source terms respectively. The VOF transport equation is given below. 
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where a = 1 stands for liquid and a = 0 is for gas. The interface is located at 1 > a > 0. S a represents the 
volume transfer rate across the two-phase boundaries. And, p m = (1— a)p g +ap, where Pm denotes the time- 

mean density of the mixture, p ( and pi denote gas and liquid density respectively, u t and u't are the i-component 
of the density-weighted mean and fluctuating part of the instantaneous velocity, 4) and are the density-weighted 
mean and fluctuating part of the instantaneous scalar quantities including the species concentrations, turbulence 
quantities and the gas mixture enthalpy, p is the mean pressure and S is used to model sources terms due to mass 
transfer, momentum transfer, species production, etc. The mixture density is used in the governing equation to 
represent the present continuum interface model. For pure-liquid or pure-gas regions, the density would recover 
liquid density or gas density respectively based on the mixture density and the interface definitions. Detailed 



expressions of the source terms can be found in refs. 1 and 2. In the present study, the turbulence correlation 
terms, u ', uV and u\ <t>\ are modeled by the two-equation turbulence closure models 14 using eddy viscosity 

formulations. For a given solution of the a field, equation (1) can be rewritten as the following form to maintain 
accurate transient from the gas-phase flowfield to the liquid-phase flowfield through the interface. 
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CONTINUUM SURFACE TENSION FORCE MODEL 

The surface tension forces in the continuum surface force (CSF) model are formulated as continuous body 
forces across the interface. These forces can be written as. 
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where a denotes the surface tension coefficient and n stands for a unit normal vector at the interface calculated 
based on the solution of a. 


INTERFACE MASS TRANSFER MODEL 

To simulate the mass transfer effects along the interface of the gas and liquid phases, the source term of Eq. 
(2) is modeled. The sign and magnitude of the VOF source term, 5 0 in Eq. (2), depend on the physical processes 
involved. For spray atomization applications, for instance, it would represent the volume stripping rate which 
can be modeled by atomization correlations. For spray coating, condensation and chemical vapor deposition 
processes, on the other hand, positive volume flow rates would result 

There are basically two ways of modeling the VOF source term. One method is to treat the inter-phase 
volume flow rate as convection process along the interface. That is. 


S_ = v. Va 

where V, denotes the volume exchange velocity. This method allows the source term and the convection terms be 
combined so that the source effects are handled implicitly. However, the accuracy of this approach depends 
entirely on the assignment of the volume exchange velocity. Vectors directly normal to the interface do not 
always guarantee good solution. 

The second approach, which is the most straight forward treatment, solve the VOF equation with the source 
term specified explicitly. This approach is general and provides good numerical accuracy. The only drawback of 
this method is the possible numerical instability related to high volume flow rates. Time step size must be small 
enough in order to solve the VOF equation with large source term. For shear coaxial spray, empirical 
correlations such as the Reitz model 15 can be used to calculate the mass stripping rate and spray atomization 
droplet size distributions. In order to provide better interfacial resolution which is essential for the success of the 
present VOF method, all test cases were simulated with a high-order upwind Total Variation Diminishing (TVD) 


Scheme. Only the convection terms were modeled by using the TVD flux limiters. The convection terms of the 
VOF transport equation can be expressed by finite difference approximation as: 
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where f and h represent first-order fluxes and TVD flux limiters respectively. The TVD flux limiters function as 
anti-diffusion terms to recover high-order accuracy. The first-order fluxes and the TVD flux limiters are given 
below. 
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where the minmod functions in the TVD flux limiters are written as: 
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The order of accuracy of this scheme is determined by the parameters y and p. Only the second-order upwind 
schemes was used in this study. That is, 
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where y = -1. The compression factor, p, is used to sharpen the interface for better resolution. 

COMBUSTION MODEL 

For reacting flow computations, a penalty function treatment 16, 17 of the reacting source terms is employed to 
calculate the source terms of the species equaitons. Finite-rate chemical kinetics for the 02/H2 1 * and 02/RP- 
1/H2 19 systems were employed for good representation of the reacting flow development Global reaction models 
were used during the ignition and flame spreading period to achieve better computational efficiency for spray 
combustion applicaitons. This approach is justified for steady-state calculations. 

VALIDATION 


LAMINAR IMPINGING JETS 

Impinging elements are common in injector designs. The prediction of the liquid sheet development has 
been a challenging task for most CFD methods. With the present VOF model, the problem of modeling the 
liquid jets impingement process is solved with very minimal computational effort In this test case, three 
impinging water jets with different impinging half angle (i.e. 30°, 45° and 60°) were considered. A grid size of 
61 x 41 x 41 was generated to simulate the three-dimensional flowfields. The injection speed was 10 m/s and 
the flowfields were assumed to be laminar. The calculation, which was a transient approach, was terminated at 
200 time steps before sheet breakup may occur. Figure 1 shows the predicted flowfields and liquid surface shapes 
for the three jets. A thin sheet of liquid was formed after impingement Figure 2 shows the liquid sheet 
thickness distribution function and comparisons with an exact solution given by ref. 20. Good agreement is 
revealed in the data comparisons. Future development effort shall include models for atomization and turbulence 
effects on jet impingement 



SHEAR COAXIAL LIQUID JET ATOMIZATION 


A cold flow case of coaxial liquid jet atomization reported by Pal, et al 21 was employed to test the current 
numerical method. Water and gaseous nitrogen were used as the liquid and gas in the study with injection 
velocity 14.3 m/s and 293 m/s (Mach number around 0.86) respectively. The operating conditions were 300-k 
temperature and one atmospheric pressure. A grid size of 61 x 31 was used for the computation. Droplet 
evaporation effects were considered, which is not important in this case. Figure 3 shows the velocity vectors and 
liquid core. Liquid jet velocity stays fairly constant due to large density difference between two fluids. The 
numerical model is stable even for large time-step size and incompressible (liquid) /compressible (gas) flow 
situation. Figure 4 presents the particle size variations and locations. Particle turbulent dispersion effect has 
been included and large particle dispersion is shown due to high turbulent intensity of the gas phase. The 
numerical prediction is reasonable compared to the observed main flow feature 21 . 

SHEAR COAXIAL INJECTOR SPRAY COMBUSTION 

A test case of LOX/H2 coaxial spray with the same operating conditions as the Case D of the Penn State uni- 
element experiment was investigated with the present VOF/particulate multi-phase model and a 7-step H2/02 
finite rate chemical kinetics. The LOX and GH2 flowrates were 0.373 Ibm/sec and 0.072 lbm/sec respectively. 
The velocity ratio at the injector outlet was 28.6 with the hydrogen gas going at a relatively high speed around 
Mach 0.6. The mean droplet size of the primary atomization under the operating conditions was about 30 
micron. A mesh with size 121 x 61 was used for the computation. A cold flow case was calculated for 1000 time 
steps to have a better gas mixture established before a hot spot of 1000 degree-K was introduced near the 
chamber wall in the middle of the chamber length. Once the flame was started it spread in both upstream and 
downstream directions. At the final stage the flame front was going very slowly upstream until it’s near the 
injector face. The total number of time steps for one case was around 5000 (nondimensional time step size of 
0.2). Results of the computations are shown in Figure 5 for the temperature and particle fields. These results 
give reasonable representation of the observed experimental flowfield. The liquid jet interface is also shown in 
Figure 5b. The calculated chamber pressure is a little higher (10 percent) than the measured value indicating 
larger flame zone was predicted. This is relating to the turbulence mixing problem. 

Another case for LOX/RP-1/H2 tripropellant coaxial injector case is also investigated. Same chamber 
geometry as the previous case, except for the injector details which was designed by the author, was used. The 
LOX, RP-1 and H2 flowrates were assumed to be 0.9659 lbm/sec, 0.1393 lbm/sec and 0.0495 lbm/sec 
respectively. This high flow rates produced a high chamber pressure around 38 ATM. As opposed to the 
LOX/H2 case, the ignition process for this case was some what tricky. A flame zone was finally established after 
many attempts of igniting the propellant More investigation is needed to study the key factor affecting the 
ignition process. Preliminary results of this case are shown in Figure 6. It is clear that a much longer stand-off 
distance was predicted for this case. The particle and VOF surface plot (Figure 6b) shows that the current 
multiple VOF model is working. 


CONCLUSIONS 

This paper summarizes the technical development and validation of a unified multiphase computational fluid 
dynamics (CFD) numerical method using volume-of-fluid (VOF) model which can be employed to analyze 
general multiphase flow problems involving free surface. A time-accurate multiphase viscous flow solution 
method has been developed and validated with benchmark test cases including laminar impinging jets and a 
coaxial jet spray atomization simulation. The gas-liquid interface mass, momentum and energy conservation 
properties are modeled by continuum surface force mechanisms. A new solution method is developed such that 
the present VOF model can be applied for all-speed flow regimes. The main achievements of the present study 
are the development and verification of the fractional volume-of-fluid cell partitioning approach into a predictor- 
corrector algorithm to deal with multiphase (gas-liquid) free surface flow problems; and the coupling of this 
unified all-speed algorithm with a non-iterative Eulerian/Lagrangian particle tracking and finite-rate chemistry 
model. The successful implementation of the VOF model in this study can lead a reliable numerical model for 



spray atomization and combustion flow analyses as well as a fundamental building block for a numerical tool 
which can be used to study complex multiphase flow phenomena which are sometimes bard to visualize 
experimentally. 
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Data— ID 
a 1 0599E+O2 

b 4 . 5497E+02 

c 8 0394E+02 

d 1 15238+03 

e 1 5018E+O3 

f 1 8508E+03 

g 2 1996E+03 

h 2 5488E+03 

i 2 8977E+03 

j 3 24G7E+03 

k 3 5957E+03 



Figure 5. LOX/H2 shear coaxial uni-element spray combustion test case. 

(a) Temperature contours and (b) VOF surface and particle plot 



(a) 


Data- ID 
a 1 0000E+O2 

b 4 284EE+02 

c 7 5692E+02 

d 1 0853E+03 

e 1 4138E+03 

f 1 7423E+03 

g 2 O707E+O3 

h 2 3992E+03 

i 2 7276E+03 

j 3 O561E+03 

k 3 384GE+03 
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Figure 6. L0X/RP-1/H2 shear coaxial uni-element spray combustion test case, 
(a) Temperature contours and (b) VOF surface and particle plot. 




